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The phase diagram of electron-doped pnictides is studied varying the temperature, electronic 
density, and isotropic quenched disorder strength by means of computational techniques applied to 
a three-orbital ( xz , yz, xy) spin-fermion model with lattice degrees of freedom. In experiments, 
chemical doping introduces disorder but in theoretical studies the relationship between electronic 
doping and the randomly located dopants, with their associated quenched disorder, is difficult to 
address. In this publication, the use of computational techniques allows us to study independently 
the effects of electronic doping, regulated by a global chemical potential, and impurity disorder at 
randomly selected sites. Surprisingly, our Monte Carlo simulations reveal that the fast reduction 
with doping of the Neel Tn and the structural Ts transition temperatures, and the concomitant 
stabilization of a robust nematic state, is primarily controlled by the magnetic dilution associated 
with the in-plane isotropic disorder introduced by Fe substitution. In the doping range studied, 
changes in the Fermi Surface produced by electron doping affect only slightly both critical tempera¬ 
tures. Results obtained varying the disorder strength indicate that the specific material dependent 
phase diagrams experimentally observed are a consequence of the variation in disorder profiles in¬ 
troduced by the different dopants. Our results are also in agreement with neutron scattering and 
scanning tunneling microscopy, unveiling a patchy network of locally magnetically ordered clusters 
with anisotropic shapes, even though the quenched disorder is locally isotropic. This study reveals 
a remarkable and unexpected degree of complexity in pnictides: the fragile tendency to nematicity 
intrinsic of translational invariant electronic systems needs to be supplemented by quenched disorder 
to stabilize the robust nematic phase experimentally found in electron-doped 122 compounds. 

PACS numbers: 74.70.Xa, 74.25.-q, 74.25.Dw, 71.10.Fd 


I. INTRODUCTION 

The mechanism that leads to high critical temperature 
superconductivity in iron-pnictides is still elusive, 

mainly because the several simultaneously active degrees 
of freedom (d.o.f.) in these materials pose a major theo¬ 
retical challenge. While magnetic mechanisms are often 
invoked to explain the d -wave superconductivity in the 
cuprates HE], the role of the orbitals is added to the 
mix in the case of the iron-based compounds. Moreover, 
the symmetry of their superconducting state is still under 
considerable debate [2]. 

The interaction among the many different d.o.f. in 
pnictides generates rich phase diagrams when varying 
temperature and doping [9]. In addition to the supercon¬ 
ducting phase, magnetic and nematic phases, accompa¬ 
nied by structural distortions, have been identified [M3]. 
To properly address this difficult problem it is necessary 
that the spin, orbital, lattice, and charge should all be 
incorporated in a treatable model that allows to monitor 
their respective roles in the properties of these materials. 
Due to the complexity of the problem most of the previ¬ 
ous theoretical studies have been performed either in the 
weak or strong coupling limits. In weak coupling, the in¬ 
teractions among the electrons are considered small and 
the physical properties are studied in momentum space 
in terms of itinerant electrons, with emphasis on par¬ 
ticular properties of their Fermi Surfaces (FS) such as 


nesting mm- On the other hand, the strong coupling 
approach is based on the experimental observation of lo¬ 
calized magnetic moments and on the fact that several 
properties of the pnictides can be reproduced via Heisen¬ 
berg models [I8lf2()j . Both approaches were successful 
in the study of the magnetic properties of the parent 
compounds, indicating that in these materials both lo¬ 
calized and itinerant magnetic moments are important. 
However, upon doping there are challenges explaining ex¬ 
perimental data in both approximations. In particular, 
when doping is achieved by chemical substitution of iron 
atoms then the associated effects of disorder must also 
be incorporated into the theoretical considerations. 

The parent compound of the 122 family, BaFe 2 As 2 , 
can be doped with electrons by replacing Fe by a transi¬ 
tion metal (TM) resulting in Ba(Fei_a;TM x ) 2 As 2 or with 
holes by replacing Ba by an alkali metal (A) leading to 
Bai_ 2 ,A a ;Fe 2 As 2 [8]. It is also possible to dope in an 
isovalent manner replacing, for example, Fe with Ru to 
obtain Ba(Fei_ a; Ru a ;) 2 As 2 [21]. Nominally, replacing Fe 
with Ru, Co, Ni, and Cu would introduce 0, 1, 2, and 3 
electrons per dopant atom. However, experiments indi¬ 
cate a difference between nominal doping x and the mea¬ 
sured doping concentration x m usually determined us¬ 
ing wavelength dispersive x-ray spectroscopy (WDS) [9]. 
This means that in some cases, electrons can get trapped 
by the doped impurities [22]. Chemical substitution in¬ 
troduces an amount of disorder that is difficult to control 
experimentally. In addition to electrons being trapped, 
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other effects such as magnetic dilution and impurity scat¬ 
tering may also occur [23 1 . 

In undoped 122 compounds the structural and the 
Neel transition temperatures, Tg and TV, are equal to 
each other. Upon electron doping both are rapidly re¬ 
duced, with Tg decreasing at an equal or slower rate 
than T/v B EH- The reduction of these temperatures 
is explained in weak coupling by a loss of FS nesting in¬ 
duced by the electronic doping and in strong coupling 
by magnetic dilution as in t-J models. However, these 
views seem to be in contradiction with several experimen¬ 
tal results. For example, in Ba(Fei_ a ,Ru a; ) 2 As 2 , which 
nominally does not introduce electronic doping and as¬ 
sociated changes in FS should not be expected, both Tg 
and TV decrease with doping and the material eventu¬ 
ally becomes superconducting |21j . In addition, doping 
with Co, Ni, and Cu is expected to introduce 1, 2, and 
3 extra electrons per doped atom. However, the experi¬ 
mentally observed reduction on TV and Tg was found to 
be primarily a function of the doping x rather than of the 
density of electrons [SUM!- Experiments, thus, indicate 
that when dopants are introduced directly on the Fe-As 
planes, as it is the case for electron-doped 122 materials, 
disorder must play an important role m [13 [2SH2Z1- Due 
to the experimental uncertainty on the actual doping con¬ 
centration and the nature of the disorder, a theoretical 
understanding of the phase diagrams under these chal¬ 
lenging circumstances is elusive. Density functional the¬ 
ory (DFT) studies indicated that in-plane-doped atoms 
would tend to trap electrons [22], while first-principles 
methods found that the interplay between on-site and 
off-site impurity potentials could induce FS distortions 
in nominally isovalent doping [23j . Moreover, a calcu¬ 
lation considering two-orbiton processes predicts a non- 
symmetric impurity potential which could be responsible 
for the observed transport anisotropies Kl¬ 
in this publication, the effects of electron doping in the 
122 pnictides will be studied numerically using a spin- 
fermion model (SFM) for the pnictides [28H50] including 
the lattice d.o.f. [31]. The SFM considers phenomenolog¬ 
ically the experimentally gathered evidence that requires 
a combination of itinerant and localized d.o.f. to properly 
address the iron-based superconductors [3] [3] [52] [55]. The 
itinerant sector mainly involves electrons in the xz, yz 1 
and xy e?-orbitals !-TTj while the localized spins represent 
the spin of the other d-orbitals [25J 25], or ' n a Landau- 
Ginzburg context it can be considered as the magnetic 
order parameter. 

The focus of this effort will be on the structural and the 
Neel transitions, and the properties of the resulting ne¬ 
matic phase that will be monitored as a function of the 
electronic and impurity densities. Earlier studies per¬ 
formed in the undoped parent compounds indicated that 
the coupling between the lattice orthorhombic distortion 
ei, associated to the elastic constant Cg6, and the spin- 
nematic order parameter Hq stabilizes the orthorhombic 
( 7 r, 0) antiferromagnetic (AFM) ground state [31] with 
Tg = T/v as in the 122 materials [9]. The small sepa¬ 


ration between Tg and T/v observed in the parent com¬ 
pounds of the 1111 family [35] was found to be regulated 
by the coupling of the lattice orthorhombic distortion to 
the orbital order parameter <Iq m- 

This is the first time that electronic doping is com¬ 
putationally studied in a system that includes mag¬ 
netic, charge, orbital, and lattice d.o.f. supplemented 
by quenched disorder. Our numerical approach involves 
Monte Carlo (MC) calculations on the localized spin and 
lattice components, combined with a fermionic diagonal- 
ization of the charge/orbital sector. In addition, twisted 
boundary conditions (TBC) and the Travelling Cluster 
Approximation (TCA) are implemented [36j in order to 
study large clusters of size 64 x 64, a record for the spin- 
fermion model. This numerical approach allows us to in¬ 
corporate the effects of random on-site and off-diagonal 
disorder and to obtain results for temperatures above 
Tg where all d.o.f. develop strong short-range fluctua¬ 
tions HU EZ], a regime difficult to reach by other many- 
body procedures. Our main conclusion is that quenched 
disorder is needed to enhance the (weak) electronic ten¬ 
dency to form a nematic phase in 122 materials. That 
a critical temperature such as T/v decreases faster with 
doping by including disorder than in the clean limit is 
natural [38] 139], but our most novel result is the concomi¬ 
tant stabilization of a nematic regime. In other words, 
T/v and Tg are affected differently by disorder. Isotropic 
disorder is sufficient to obtain these results. Our analy¬ 
sis illustrates the interdependence of the many degrees of 
freedom present in real materials and the need to study 
models with robust many-body techniques to unveil the 
physics that emerges in these complex systems. 

The organization of the paper is as follows: the model 
is described in Section [H and the computational methods 
are presented in Section Hi] Section [TV] is devoted to the 
main results addressing the phase diagram upon doping. 
Section [V] describes the properties of the nematic phase 
stabilized in our study, including a comparison with neu¬ 
tron scattering and scanning tunneling microscopy ex¬ 
periments. The discussion and summary are the scope of 
Section [VI] 


II. MODEL 
A. Hamiltonian 

The spin-fermion model Hamiltonian studied here is 
based on the original purely electronic model [25M50] sup¬ 
plemented by the recent addition of couplings to the lat¬ 
tice degrees of freedom [55]: 

FIsF — ^Hopp T -^Hurid T -^Heis “b ^SL T ■ (1) 

Haopp is the three-orbitals ( d xz , d yz , d xy ) tight-binding 
Fe-Fe hopping of electrons, with the hopping ampli¬ 
tudes selected to reproduce ARPES experiments. Read¬ 
ers can find these amplitudes in previous publications, 
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such as in Eqs.(l-3) and Table 1 of Ref. [33]. The av¬ 
erage density of electrons per iron and per orbital is 
n= 4/3 in the undoped limit [333 and its value in the 
doped case is controlled via a chemical potential in¬ 
cluded in -f/nopp [SB]. The Hund interaction is stan¬ 
dard: i?Hund=—J h Si a ' s i,a) with Si the localized 
spin at site i and s, a the itinerant spin corresponding 
to orbital a at the same site [401. Hn els contains the 
Heisenberg interaction among the localized spins involv¬ 
ing both nearest-neighbors (NN) and next-NN (NNN) in¬ 
teractions with respective couplings Jnn and Jnnn> and 
a ratio Jnnn/>/nn = 2/3 (any ratio larger than 1/2 would 
have been equally effective to favor “striped” spin order). 
Having NN and NNN Heisenberg interactions of compa¬ 
rable magnitude arise from having comparable NN and 
NNN hoppings, caused by the geometry of the material. 

The coupling between the spin and lattice degrees of 
freedom is given by HsL=—g ’Jqei -TS) :TB], where g is 
the spin-lattice coupling [41]. Finally, Hsus is the spin 
stiffness given by a Lennard-Jones potential that speeds 
up convergence, as previously discussed [3B]. Note that 
the lattice-orbital coupling term, Hol =—A JA ^iCi j3B] . 
is omitted because previous work indicated that A in¬ 
duces a (small) nematic phase with Tg > Tjy directly 
in the parent compounds ED EB]. Since the goal of the 
present effort is to study the 122 family, characterized 
by T s = Tjv in the undoped case, then this term is not 
included to reduce the number of parameters. 



FIG. 1: (color online) Internal structure of dopant sites. 
Sketch shows the location of a dopant where the magnitude 
of the localized spin, Si, is reduced from the original value S. 
In addition, the neighboring localized spins are also assumed 
to be affected by the presence of the dopant. The four imme¬ 
diate nearest-neighbors have a new localized spin magnitude 
Snn, while the four next nearest-neighbors have a new local¬ 
ized spin magnitude Snnn, such that Si < Snn < Snnn < S 
(S is the undoped localized spin magnitude, assumed to be 1 
in this publication unless otherwise stated). 


B. Quenched Disorder 

On-site diagonal disorder is introduced by adding an 
impurity potential //(id) to Ni randomly selected sites 
id where transition metal atoms replace Fe. The density 
of impurity atoms x is defined as x = Ni/N, where N is 
the total number of lattice sites. In addition, the value 
of the localized spin at the impurity site, Si, is reduced 
since, for example, Co dopants in BaFe 2 As 2 are non¬ 
magnetic [42]. This effectively reduces the local Hund 
coupling Jhj and the spin-lattice coupling g\ at the im¬ 
purity sites. We also will study the effect of extending 
the spatial range of the impurity by reducing the values 
of the localized spins to Snn (Snnn) at the NN (NNN) 
of the impurity sites with the corresponding effective de¬ 
crease in Jh and g at those sites (see Fig. [T]). Thus, 
off-diagonal isotropic disorder results from the effective 
reduction of the Heisenberg couplings at the bonds con¬ 
necting the impurity sites and their neighbors [40]. Note 
that off-diagonal disorder could also be introduced in the 
eight hopping amplitudes present in Hn opp [361 but for 
simplicity we decided not to consider hopping disorder at 
this time. 


III. METHODS 

The Hamiltonian in Eq.Q was studied via a well- 
known Monte Carlo method [30l 03j applied to (i) the 
localized (assumed classical) spin degrees of freedom Si 
and (ii) the atomic displacements that determine the lo¬ 
cal orthorhombic lattice distortion ei mm- For each 
Monte Carlo configuration of spins and atomic positions 
the remaining quantum fermionic Hamiltonian is diago¬ 
nalized. The simulations are performed at various tem¬ 
peratures, dopings, and disorder configurations and lo¬ 
cal and long-range observables are measured. Note that 
with the exact diagonalization technique results can be 
obtained comfortably only on up to 8 x 8 lattices, which 
may be too small to provide meaningful data at the low 
rates of doping relevant in the pnictides. For this rea¬ 
son we have also used the Traveling Cluster Approxima¬ 
tion 03] where a larger lattice (64x64 sites in most of this 
effort) can be studied by performing the MC updates via 
a travelling cluster centered at consecutive sites i, with a 
size substantially smaller than the full lattice size of the 
entire system. Twisted boundary conditions were also 
used [45] to obtain (almost) a continuum range of mo¬ 
menta. For simplicity, most couplings are fixed to values 
used successfully before [30]: Jh=0.1 eV, Jnn=0.012 eV, 
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and Jnnn=0.008 eV. The dimensionless version of the 
spin-lattice coupling g is fixed to 0.16 as in [3l|. The fo¬ 
cus of the publication is on the values for the parameters 
associated with disorder and the corresponding physical 
results, as discussed in the sections below. 

An important technical detail is that to improve nu¬ 
merical convergence, and to better mimic real materials 
that often display an easy-axis direction for spin orien¬ 
tation, we have introduced a small anisotropy in the x 
component of the super-exchange interaction so that the 
actual Heisenberg interaction is: 

Hue is = Jnn 5> • s i + SS i S D 

" j> ( 2 ) 

+Jnnn ^ (Si-S m + SSfS^), 

«im» 

with S = 0.1. This anisotropy slightly raises T/v, but 
the magnetic susceptibility X S becomes much sharper at 
the transition temperatures, facilitating an accurate de¬ 
termination of T/v. 

The Monte Carlo simulations with the TCA procedure 
were mainly performed using 64 x 64 square lattices |46| . 
Typically 5,000 MC steps were devoted to thermaliza- 
tion and 10,000 to 25,000 steps for measurements at each 
temperature, doping, and disorder configuration. The 
results presented below arise from averages over five dif¬ 
ferent disorder configurations. The expectation values 
of observables remain stable upon the addition of extra 
configurations due to self-averaging. The magnetic tran¬ 
sition was determined by the behavior of the magnetic 
susceptibility defined as 

Xs(k,o) = N(3(S(tt,0) — (5(7r, 0))) 2 , (3) 

where j3 = 1/fcgT, N is the number of lattice sites, and 
S(-k,0) is the magnetic structure factor at wavevector 
(7r, 0) obtained via the Fourier transform of the real-space 
spin-spin correlations measured in the MC simulations. 
The structural transition is determined by the behavior 
of the lattice susceptibility defined by 

Xs =Np(6-(d)) 2 , (4) 

where (5 = (fAzlhi) an d ai j s the lattice constant along 
the i = x or y directions. These lattice constants are 
determined from the orthorhombic displacements e\ [56] . 


IV. RESULTS 

Our first task is to understand the effect of doping and 
disorder on the magnetic and structural transitions. For 
this purpose, we studied the evolution of T)v and Tg vs. 
doping concentration under different disorder setups. 


A. Clean limit 


Consider first the “clean limit”. The red squares in 
Fig- E show the evolution of Tiv and Tg when the elec¬ 
tronic doping does not introduce disorder. In this case 
is hardly affected and it continues to be equal to Tg 
for all dopings investigated here. This result indicates 
that the reduction of T/v and Tg, and the stabilization 
of a nematic phase in between the two transitions ob¬ 
served experimentally upon electron doping [5] , does not 
emerge just from the reduction of Fermi Surface nest¬ 
ing induced by the electronic doping. This conclusion is 
not surprising if we recall that the undoped V-site lat¬ 
tice has 4 N electrons which means that for x = 10% 
the number of added electrons is N e = 0.1 N and, thus, 
the percentual change in the electronic density is just 
100 x (0.1N/4N) = 2.5%. Such a small percentual varia¬ 
tion in the electronic density should not produce substan¬ 
tial modifications in the FS, explaining why the changes 
in nesting are small and, thus, why the critical tempera¬ 
tures are not significantly affected. Then, disorder maybe 
needed to understand the experiments. 



FIG. 2: (color online) Clean limit and effect of Co doping. 
The clean limit results (open and solid red points) indicate 
that T s = Tn and both are approximately constant in the 
range studied. For Co doping, the Neel temperature Tv (open 
circles and black dashed line) and the structural transition 
temperature Ts (filled circles and black solid line) vs. the 
percentage of impurities x are shown. The on-site disorder 
is I\ = —0.1 and the off-diagonal disorder is determined by 
Si = 0, Snn = S/4, and Snnn = S/2. For both sets of curves, 
i.e. with and without quenched disorder, the density of doped 
electrons equals x to simulate Co doping. The cluster used 
has a size 64 x 64. 














5 


B. Co doping 

To study the effect of quenched disorder, let us first fo¬ 
cus on Co doping, which nominally introduces one extra 
electron per dopant. In Fig. [2j the Neel and structural 
transition temperatures are presented for the case where 
one extra electron is contributed by each replaced iron 
atom, which means that x = n, where n is the density 
of added electrons and x is the density of replaced iron 
atoms. We considered several possible values for the on¬ 
site impurity potential and spin values near the impurity 
(see details discussed below) and we found that the exper¬ 
imental data of Ref. H]were best reproduced by setting the 
on-site impurity potential as Ij = — 0.1 (in eV units) [47] 
and by using Si = 0 at the impurities since there is evi¬ 
dence that Co doped in BaFe 2 As 2 is non-magnetic 02] , 
This effectively sets to zero the Hund coupling Jh,i and 
the spin-lattice coupling g\ at the impurity sites. In ad¬ 
dition, we also reduced the localized spins to 5/4 (5/2) 
at the NN (NNN) of the impurity sites with the corre¬ 
sponding effective decreased in Jp and g at those sites. 
The overall chemical potential g was adjusted so that the 
density of added impurities equals the density of added 
electrons, which corresponds to an ideal Co doping [5]. 



T (K) 

FIG. 3: (color online) The magnetic susceptibility (open black 
symbols) and the lattice susceptibility (filled red symbols) vs. 
temperature. The sharp peaks indicate the Neel temperature 
Tjv and the structural transition temperature Ts for the case 
of 5% Co-doping. The on-site disorder is I\ = —0.1 and the 
off-diagonal disorder is defined by Si = 0, Snn = S/4, and 
Snnn = S/2. The cluster used is 64 x 64. 


The black filled (open) circles in Fig.[2]show the evolu¬ 
tion with impurity doping of the structural (Neel) transi¬ 
tion temperatures in the presence of the disorder caused 
by replacing Fe by Co at random sites. The magnetic 
dilution induced by doping causes a rapid reduction in 
Ts and Tjv, similarly as observed in experiments [9], 


and remarkably also opens a robust nematic phase for 
Tjv < T < Ts since disorder affects differently both tran¬ 
sition temperatures. In fact, the separation between T n 
and Ts is very clear in the magnetic and lattice suscepti¬ 
bilities that are displayed for 5% doping, as example, in 
Fig. m The magnetic properties of the different phases 
are also clear by monitoring the behavior of the real- 
space spin-spin correlation functions presented in Fig. [4] 
In panel (a) for T = 120 K (T > Ts) the spin correlations 
effectively vanish at distances larger than two lattice con¬ 
stants and there is no difference between the results along 
the x and y axes directions, indicating a paramagnetic 
ground state. However, at T = 95 K (Tjv < T < Ts), 
panel (b), the correlations now display short-range AFM 
(FM) order along the x ( y) directions demonstrating the 
breakdown of the rotational invariance that characterizes 
the nematic phase, but without developing long-range or¬ 
der as expected. Lowering the temperature to T = 80 K 
(T < Tjv), panel (c), now the correlations have devel¬ 
oped long range (7r,0) order, as expected in the antifer¬ 
romagnetic ground state. To our knowledge, the results 
in figures such as Fig.[2]provide the largest separation be¬ 
tween Ts and Tjv ever reported in numerical simulations 
of realistic models for iron-based superconductors. 



FIG. 4: (color online) Real-space spin-spin correlation func¬ 
tions vs. distance on a 64 x 64 lattice; (a) corresponds to 
T = 120 K (T > Ts) in the paramagnetic regime, (b) to 
T = 95 K (Tjv < T < Ts) in the nematic state, and (c) to 
T = 80 K (T < Tv) in the long-range ordered magnetic state. 
The AFM correlations along x are indicated with solid cir¬ 
cles while the FM correlations along y are denoted with open 
circles. The results are for 5% Co-doping with off-diagonal 
disorder set by Si = 0, 5nn = 5/4, and Snnn = 5/2. 
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C. Cu doping 

Let us consider now the effect of doping with Cu which, 
nominally, introduces three electrons per doped impu¬ 
rity For this purpose we increased the chemical po¬ 
tential at a faster rate so that the added density of elec¬ 
trons is n = 3.t, instead of n = x as for Co doping. The 
results are shown in Fig. [5] When the critical tempera¬ 
tures for both Cu and Co doping are plotted as a function 
of the density of impurities x, in Fig. [5ja) it can be seen 
that the results are approximately independent of the 
kind of dopant. This indicates that the critical tempera¬ 
tures are primarily controlled by the amount of quenched 
disorder (namely, by the number of impurity sites) rather 
than by the actual overall electronic density, at least in 
the range studied. This conclusion is in excellent agree¬ 
ment with the experimental phase diagrams shown, for 
example, in Fig. 26(a) of Ref. [3], for the case of several 
transition metal oxide dopants. Thus, working at a fixed 
electronic density n, the values of Ty and T$ are smaller 
for the case of Co doping than for the case of Cu-doping, 
as shown in Fig. [5jb) , because more Co than Cu impu¬ 
rities have to be added to achieve the same electronic 
density, underlying the fact that Co doping introduces 
more disorder than Cu doping at fixed n. These results 
are also in good agreement with the experimental phase 
diagram in Fig. 26(b) of Ref. [9). 

D. Dependence on impurity characteristics 

Let us consider the dependence of the Neel and the 
structural transitions temperatures on the local details of 
the magnetic dilution caused by the disorder. In Fig. [6] 
results for Tjv and T$ are shown as a function of impurity 
doping with the chemical potential set to introduce one 
electron per dopant. The clean limit data (red squares, 
case I) is displayed again for the sake of comparison. The 
blue triangles (case II) are results for I\ = —0.1 and 
Si = S/2, leaving Snn and Snnn untouched (i.e. equal 
to S). This ultra local magnetic dilution induces effec¬ 
tive NN and NNN reductions in the Heisenberg couplings 
accelerating the rate of decrease of the critical tempera¬ 
tures. However, the nematic phase is still not stabilized 
and, thus, it does not reproduce the experimental be¬ 
havior for the Co-doped parent compound. Reducing Si 
to zero, as indicated by the green diamonds in the figure 
(case III) and keeping Snn and Snnn untouched, slightly 
increases the rate of reduction of the critical tempera¬ 
tures with doping and stabilizes the nematic phase only 
after a finite amount of doping x ~ 10% has been added 
but in a very narrow range of temperature. The conclu¬ 
sion of cases I, II, and III is that a very local description 
of the dopant is insufficient to reproduce experiments. 

We have found that in order to generate a robust ne¬ 
matic phase upon doping, extended effects of magnetic 
dilution must be considered. The upside-down purple 
triangles (case IV) in Fig. [6] show results for Si = S/2, 




FIG. 5: (color online) Contrast of effects of Cu and Co doping. 
The Neel temperatures T,v (dashed lines) and the structural 
transition temperatures Ts (solid lines) for Co doping (black 
open and solid circles) and for Cu doping (blue open and solid 
triangles) are shown. Results are presented first (a) vs. the 
impurity density x and second (b) vs. the added electronic 
density n. The off-diagonal disorder is set at Si = 0, Snn = 
S/4, and Snnn = S/2. The cluster size is 64 x 64. 


Snn = 0.7S, and Snnn = 0.9S. The nematic regime 
is still too narrow. But the results for Si = 0 with 
Snn = S/4 and Snnn = S/2 (black circles, case V), 
already shown in Fig. [2j indicate that increasing the 
strength of the extended off-diagonal disorder does in¬ 
duce a faster reduction of the critical temperatures and 
stabilizes a larger nematic region. Our computer simu¬ 
lations suggest that the range and strength of disorder, 
specifically the extended magnetic dilution, is crucial for 
the stabilization of the nematic phase when Tty = T<j in 
the parent compound. 
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We have observed that the effect of the on-site impu¬ 
rity potential I\ is weak. In principle, we could have 
kept the overall chemical potential /i fixed and control 
the added electronic density n by merely adjusting the 
values of the impurity potential. However, this does not 
induce noticeable changes in the critical temperatures, 
due to the small overall modifications in the electronic 
density discussed before. This is not the manner in which 
doping seems to act in the real electron-doped pnictides. 
Thus, we believe that working with a fixed value of the 
impurity potential and adjusting the electronic density 
with the overall chemical potential allows to study the 
effects of isotropic quenched disorder and varying elec¬ 
tronic density in a more controlled and independent way. 

Considering the negligible effect on the critical tem¬ 
peratures caused by pure electronic doping (clean limit) 
and, by extension, the on-site impurity potential, the re¬ 
sults in Fig.[6]shed light on the case of isovalent doping in 
which Fe is replaced by Ru. This procedure introduces 
disorder but, at least nominally, no electronic doping. 
Experimental efforts have observed that in this case Tjv 
and Tg still decrease with doping, despite no apparent 
changes in the Fermi surface, but at a slower rate than 
with non-isovalent doping. Moreover, the critical tem¬ 
peratures do not separate from each other, i.e., no ne¬ 
matic phase is stabilized [51]. Our results lend support 
to the view that the decrease of T/v and Tg observed 
with Ru-doping is mainly due to the magnetic dilution 
introduced by doping rather than by more subtle effects 
on the electronic density which in turn would affect the 
nesting of the FS [21, . Experiments have determined 

that doped Ru is magnetic [35] which would translate to 
larger values of Si, SnNi and S'nnn in our model. In fact, 
the blue triangles (case II) in Fig. [^qualitatively capture 
the slower decrease rate and negligible separation with 
impurity doping for T/v and Tg experimentally observed 
for Ru doping [21]. 

V. PROPERTIES OF THE NEMATIC PHASE 

Having stabilized a robust nematic regime, let us study 
its properties. 


A. Neutron scattering 

Considering the importance of neutron scattering ex¬ 
periments in iron superconductors, we studied the elec¬ 
tronic doping dependence of the magnetic structure fac¬ 
tor S (k) obtained from the Fourier transform of the real- 
space spin-spin correlation functions displayed in Fig. [3] 
Experiments indicate that the low-temperature magnetic 
phase below Tg = T/v in the parent compound develops 
long range AFM (FM) order along the long (short) lattice 
constant direction in the orthorhombic lattice. This re¬ 
sults in a sharp peak at k = (n, 0) (or at (0, n) depending 
on the direction of the AFM order) that forms above the 



FIG. 6: (color online) Dependence of results with impurity 
characteristics. The Neel transition temperature Tn (dashed 
lines) and the structural transition temperature Ts (solid 
lines) vs. the percentage of impurities x for different settings 
of the off-diagonal disorder. Case I corresponds to the clean 
limit with no impurity sites (red squares). Case II has Si=S/2 
and Snn=Snnn=S untouched (blue triangles). This case may 
be sufficient for Ru doping, which is magnetic. Case III has 
Si=0 and S'nn=<S'nnn=S' untouched (green diamonds). Case 
IV has Si=S/2, Snn=0-7S, and 5nnn=0.9S i (purple upside- 
down triangles). Finally, Case V has <Si=0, iSnn=<S , /4, and 
Snnn=£'/2 (black circles). Case V appears to be the best to 
describe experiments for non-magnetic doping. The density 
of doped electrons equals x as in Co doping. In all cases the 
on-site disorder potential is kept fixed at Ji = —0.1. The 
lattice size is 64 x 64. 


small spin-gap energy [5j. More importantly for our dis¬ 
cussion and results, upon electron-doping the (7r, 0) neu¬ 
tron peak becomes broader along the direction transver¬ 
sal to the AFM order in the whole energy range [3], cre¬ 
ating an intriguing transverselly elongated ellipse. 

The results obtained numerically for 5% Co-doping 
are shown in Fig. [7] for T = 120 K (T > Tg ), i.e. in 
the paramagnetic phase. In panel (a) peaks in the spin 
structure factor S'(k) (that represents the integral over 
the whole energy range of the neutron scattering results) 
with similar intensity at wavevectors ( 7 r, 0) and (0,7r) can 
be observed. Both of these peaks are elongated along 
the direction transversal to the corresponding spin stag¬ 
gered direction, in agreement with neutron scattering [8}. 
Our explanation for these results within our spin-fermion 
model is not associated with Fermi Surface modifications 
due to electron doping, since the percentual doping is 
small as already discussed, but instead to the develop¬ 
ment of spin-nematic clusters, anchored by the magneti¬ 
cally depleted regions that form at the impurity sites. A 
Monte Carlo snapshot of the spin-nematic order parame¬ 
ter 'Iq on a 64 x 64 lattice is shown in panel (b) of Fig. [ 7 ] 
Since T > Tg, patches with (n, 0) and (0,7r) nematic or¬ 
der, indicated with green and orange in the figure, coexist 
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FIG. 7: (color online) Magnetic and nematic order in the 
paramagnetic regime. The results are for 5% Co-doping at 
T = 120 K (T > Ts) and using a 64 x 64 lattice, (a) The 
magnetic structure factor S(k), showing that the wavevectors 
(n, 0) and (0, n) have similar intensity, (b) Monte Carlo snap¬ 
shot of the spin-nematic order parameter with approximately 
the same amount of positive (green) and negative (orange) 
clusters. The impurity sites are indicated by black dots. 


in equal proportion. By eye inspection, we believe that 
the ( 7 r, 0) patches tend to be slightly elongated along the 
x direction, while the (0,7r) patches are elongated along 
the y direction. This asymmetry could be the reason for 
the shape of the peaks in the structure factor displayed 
in panel (a), since elliptical peaks can be associated to 
different correlation lengths along the x and y axes. In 
Fig . [7](a) the elliptical (-7T, 0) peak has a correlation length 
larger along the x axis than the y axis. 

The results corresponding to lowering the tempera¬ 
ture into the nematic phase (T = 95 K) are presented 
in Fig. [8] In this case the subtle effects already ob¬ 
served in the paramagnetic phase are magnified. In panel 
(a), it is now clear that the peak at (n, 0) has devel¬ 
oped a much larger weight than the peak at (0,7r), as 
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FIG. 8: (color online) Magnetic and nematic order in the ne¬ 
matic regime. The results are for 5% Co-doping at T = 95 K 
(Tjv < T < Ts) and using a 64 x 64 lattice, (a) The magnetic 
structure factor S{ k) is shown, with clear dominance of the 
(7r, 0) state, (b) Monte Carlo snapshot of the spin-nematic or¬ 
der parameter. Impurity sites are indicated by black dots. A 
positive nematic order (green) dominates, but there are still 
small pockets of negative order (orange), (c) Monte Carlo 
snapshot displaying the on-site component along the easy 
axis, S e , of the localized spin multiplied by the factor (—l) 1 ", 
with i x the axaxis component of the location of site i. Both 
the dominant blue and red clusters indicate regions with local 
( 7 r, 0) order, but shifted by one lattice spacing. This shift sup¬ 
presses long-range order when averaged over the whole lattice, 
but short-range order remains. Impurity sites are denoted as 
black dots. 
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expected. Moreover, the elongation along the transver¬ 
sal direction already perceived in the paramagnetic state 
is now better developed. The Monte Carlo snapshot of 
the spin-nematic order parameter in panel (b) shows that 
the (7r,0) (green) regions prevail over the (0,7r) (orange) 
regions, indicating that the symmetry under lattice ro¬ 
tations in the nematic phase is spontaneously broken. 
In addition, now the elongated shape of the (ir, 0) green 
clusters along the AFM direction is more clear to the eye. 
But despite the prevalence of (7T,0) clusters the system 
does not develop long-range magnetic order (compatible 
with panel (b) of Fig. [4]). This is because the many ( 7 r, 0) 
clusters are actually “out of phase” with each other. This 
is understood via the visual investigation of Monte Carlo 
snapshots, as in panel (c) of Fig. |8j where it is shown 
the component of the localized spins along the easy axis, 
S e , multiplied by a factor (—l) lx (see definition in cap¬ 
tion; the location of the impurities is indicated with black 
dots). The abundant red and blue patches all indicate 
clusters with local ( 7 r, 0) nematic order, but shifted one 
with respect to the other by one lattice spacing. The very 
small regions with (0, n) order, as in the orange regions 
of panel (b), can be barely distinguished in panel (c) with 
a checkerboard red/blue structure. 


B. Scanning Tunneling Microscopy 

The real space structure of the ( 77 , 0) nematic clus¬ 
ters obtained numerically, with an elongation along the 
x axis, can be contrasted with Scanning Tunneling Mi¬ 
croscopy (STM) measurements. In fact, STM studies 
of Co-doped CaFe 2 As 2 at 6% doping 26] [27] have al¬ 
ready revealed the existence of unidirectional electronic 
nanoestructures. These STM structures appear to have 
an average length of about eight lattice spacings along 
the AFM direction and it was argued that they may be 
possibly pinned by the Co atoms. The picture of elon¬ 
gated structures along the x axis is consistent with our 
results, as shown in panel (b) of Fig. [8j However, in our 
simulation the nematic structures are mainly located in 
between, rather than on top, the Co dopants. In our 
case this arises from the fact that the effect of disorder 
considered here reduces drastically the magnetic interac¬ 
tions at the Co or Cu dopant sites because they are not 
magnetic. 

A recently discussed new perspective is that the ne¬ 
matic state could be a consequence of anisotropic dopant- 
induced scattering rather than an intrinsic nematic elec¬ 
tronic state [2E 32], by studying the anisotropy in the 
optical spectrum [25] and in the in-plane resistivity [49] 
varying Co doping in BaFe 2 As 2 - The main argument 
to attribute the observed anisotropies to extrinsic effects 
of Co doping is that the anisotropy increases with dop¬ 
ing despite the fact that the spin order weakens and the 
lattice orthorhombicity diminishes. Our results, by con¬ 
struction, were obtained with impurity profiles that are 
symmetric under rotations of the lattice, so nematicity 


is not induced by asymmetric Co doping characteristics. 
However, we agree with the above described experimen¬ 
tal observations that quenched disorder introduced by 
the dopants is crucial for the stabilization of the nematic 
phase, otherwise in the “clean limit” there is no difference 
between Tg and Ty as already explained. 

In our simulation, the nematic phase develops be¬ 
cause the in-plane dopants allowed the formation of cigar¬ 
shaped nematic domains. These domains have shifts in 
their respective AFM orders, as it can be seen in panel 
(c) of Fig. [8j For the 122 compounds, the dopants en¬ 
hance the (weak) electronic tendency to nematicity, while 
according to our previous calculations [3T] in the par¬ 
ent compound of materials in the 1111 family, such as 
ReFeAsO (Re= La, Nd, Srn), a small temperature range 
of nematicity can be provided by the coupling between 
the lattice and the orbital degrees of freedom. This view 
may be supported by studies of the phonon modes in 
the 1111 materials m- Note also that atomic-resolution 
variable-temperature Scanning Tunnelling Spectroscopy 
experiments performed in NaFeAs, which has Tg > T)v, 
and in LiFeAs, which does not develop neither magnetic 
order nor a structural transition, indicate that cigar-like 
nematic domains develop in the nematic phase of NaFeAs 
regardless of the symmetry of the impurities observed in 
the samples m- 


VI. DISCUSSION AND CONCLUSIONS 

In this publication, the effects of electron doping in 
materials of the 122 family, such as BaFe 2 As 2 , have been 
investigated via numerical studies of the spin-fermion 
model, including charge, orbital, magnetic, and lattice 
degrees of freedom. These materials are electron doped 
via the in-plane replacement of iron atoms by transi¬ 
tion metal oxides, introducing disorder in the iron lay¬ 
ers. The results of our study suggest that the experimen¬ 
tally observed reduction of the magnetic and structural 
transition temperatures upon doping, in such a manner 
that Tpf < Tg, is primarily triggered by the influence of 
quenched disorder associated with the chemical substitu¬ 
tion of magnetic Fe atoms by non-magnetic dopants such 
as Co 22 and Cu 52] . More specifically, reducing the 
magnitude of the localized spins at and near the dopants 
rapidly reduces the values of both transition critical tem¬ 
peratures. A “patchy” nematic phase is stabilized, which 
is characterized by a majority of clusters with (7r, 0) or¬ 
der. These patches have out-of-phase magnetic order 
separated by non-magnetic regions anchored by the im¬ 
purities. While the tendency to nematicity is already 
a property of the purely electronic spin-fermion model, 
as already discussed in previous studies [31], it seems 
that for the 122 materials this fragile tendency would 
not materialize into a robust nematic phase without the 
influence of disorder. Compatible with this conclusion, 
BaFe 2 (Asi_ a: P ;c )2 (considered among the “cleanest” of 
doped pnictides since, for example, quantum oscillations 
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were observed J53]) has a splitting between Tg and TV 
which is very small (if any). 

Note that a mere change in chemical potential to in¬ 
crease the electronic doping, without adding quenched 
disorder, does not stabilize a nematic regime and intro¬ 
duces at best a very small decrease in the transition tem¬ 
peratures. This indicates that nesting effects do not play 
a major role in the opening of a robust nematic window 
with doping in 122 materials. Our results also explain the 
slower decrease of the critical temperatures, and lack of 
separation between TV and Tg, observed upon Ru dop¬ 
ing. In this case experiments have shown that Ru dopants 
in 122 materials are magnetic [55], contrary to the non¬ 
magnetic nature of Co and Cu dopants. Thus, in our 
study the values of the Hund and Heisenberg couplings 
would have to be only slightly reduced at the impurity 
sites. As shown in Fig. [6j this will reduce the rate of 
decrease, as well as the separation, of Tjv and Tg. The 
same effect may explain why TV = Tg and the decrease 
rate is slower in hole doped systems where the holes are 
introduced by replacing Ba atoms reducing the effects of 
disorder directly in the iron layers. 

In addition, the observed clusters are elongated along 
the AFM direction in agreement with similar observa¬ 
tions in STM experiments. Within the spin-fermion 
model, the cigar-like shape of the clusters is justified 
because the nearest-neighbor couplings are AFM and, 
thus, fluctuations are expected to be larger along the 
FM (frustrated) direction which reduces the associated 
correlation length. Another consequence of this behavior 
is the oval shape observed for the weight distribution of 


the magnetic structure factor around the momenta ( 7 r, 0) 
and (0,7r) for T > Tjv, in agreement with the distribution 
observed in the electron-doped case in neutron scattering 
experiments. 

In summary, we report the first computational study of 
a realistic model for pnictides that reproduces the rapid 
drop of T/v and Tg with the chemical replacement of Fe 
by transition metal elements such as Co or Cu. Since 
disorder affects differently TV and Tg, a nematic regime 
is stabilized. The key ingredient is the introduction of 
quenched disorder affecting several neighbors around the 
location of the dopant. Fermi Surface nesting effects were 
found to be too small to be the main responsible for the 
fast drop of critical temperatures. Our results are in 
agreement with neutron scattering and also with Scan¬ 
ning Tunneling Microscopy that unveiled the presence 
of anisotropic nanoclusters associated with the nematic 
state. Considering the present results for doped systems, 
together with the previously reported results for the par¬ 
ent compounds, we conclude that the spin-fermion model 
captures the essence of the magnetic properties of the 
pnictide iron superconductors. 
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